Andreas Novotny, Baptiste Serandour, Susanne Kortsch, Benoit Gauzens, Kinlan Mehdi Goulwen Jan, and Monika Winder

Contact: Andreas Novotny,

library(tidyverse)
library(xts)
library(igraph)
library(svglite)
library(circlize)
library(lubridate)
library(forcats)
library(fluxweb)
library(car)

1 Data Preparations

Figure: Overview of data preparation methods

1.1 Selectivity index based on DNA metabarcoding

Individuals of zooplankton were sampled at three locations in the Baltic Sea (Landsort Deep, Gotland Deep, and Bornholm Deep) over the season between 2017 and 2020 (Fig. E1, Tab. E3). Zooplankton samples were collected with vertical hauls from 0-30 and 30-60 m using a 90 µm-WP2 plankton net. For estimating available prey composition, water samples were collected with 10L Niskin bottles with 5 m depth intervals above the thermocline of 30 m depth or with a 20 m long hose (HELCOM 2017a). The depths were mixed by adding an equal volume of water from the Niskin bottles. 1-3L were sequentially filtered onto 25 mm diameter filters with 20 µm, 2 µm (polycarbonate), and 0.2 µm (nylon) pore sizes. Sampling for DNA analysis followed the protocol of the Swedish national monitoring program for microscopic count data.

We used a metabarcoding assay of the V3-V4 region of 16S rRNA gene of DNA extracted from the dominant zooplankton species, including Acartia spp., Centropages hamatus, Evadne nordmanni, Pseudocalanus spp., Temora longicornis, Synchaeta baltica, Bosmina spp., Keratella spp., and the water samples. Libraries were constructed using primers 341F-805R (Herlemann et al. 2011) and the sequences were assigned to the SILVA (Balvočiūtė and Huson 2017)and PhytoREF databases (Decelle et al. 2015).

Load DNA dataset zooplankton-consumer and prey availability data:

zoopl_16S <- readRDS("./Data/16SrRNA_zp.rds") %>% 
  dplyr::mutate(StationDate = paste(STATION_ID, SAMPLE_date),
                StationDateOTU = paste(STATION_ID, SAMPLE_date, OTU),
                NODE_NAME = ifelse(SORTED_genus == "S.baltica", "Synchaeta", SORTED_genus))

head(zoopl_16S)

Create table of samples, to get an overview of the number of sample replicates per day, station, and zooplankton species.

zoopl_16S %>%
  filter(SORTED_genus != "Water") %>% 
  group_by(Sample, SORTED_genus, SAMPLE_date, STATION_ID) %>% 
  summarize() %>%
  group_by(SORTED_genus, SAMPLE_date, STATION_ID) %>% 
  summarize(n = length(Sample)) %>% 
  write_csv("./Outputs/Tab.S3.csv")
## `summarise()` has grouped output by 'Sample', 'SORTED_genus', 'SAMPLE_date'.
## You can override using the `.groups` argument.
## `summarise()` has grouped output by 'SORTED_genus', 'SAMPLE_date'. You can
## override using the `.groups` argument.

To assess the selectivity index S(i) for each prey taxon (i) of a consumer, we use the standardized forage ratio as described by (Chesson 1983). S(i) = Rg(i)/Re(i)/sum:k(Rg(k)/Re(k)) Rg(i) is the Read count of i (Abundance) in the consumers gut . Re(i) is the Read count of i in the environmental sample.

Create a paired data set, with each OTU in both the water and gut sample on one row.

environment <- zoopl_16S %>% 
  # Select Water samples only
  dplyr::filter(NODE_NAME == "Water")


paired_data <- zoopl_16S %>% 
  
  # Select gut samples only
  dplyr::filter(NODE_NAME != "Water") %>% 
  
  # Merge the Zooplankton and the water dataset, on the OTU level.
  dplyr::left_join(environment, by="OTU") %>% 
  
  # Remove OTUs with no read in the water sample.
  dplyr::filter(is.na(Sample.y)==FALSE,
         Abundance.y>0)

head(paired_data)

Calculate Selectivity indices for each OTU and sample

selectivity_data <- paired_data %>% 
  
  # Calculate preference value based on water Rg(i)/Re(i)
  dplyr::mutate(Preference = Abundance.x/Abundance.y) %>% 
  
  # Summarize per sample Preference/sum(prefeerence) ; RA(ij)/RA(iw) / sum(k, RA(kj)/RA(kw))
  dplyr::group_by(Sample.x) %>%    
  dplyr::summarize(Selectivity_index=Preference/sum(Preference),
            OTU = OTU, NODE_NAME = NODE_NAME.x,
            Domain=Domain.x, Supergroup=Supergroup.x, Class = Class.x,
            Order=Order.x, Family=Family.x, Genus=Genus.x, Species=Species.x, SAMPLE_date=SAMPLE_date.x)

head(selectivity_data)

Summarize Selectivity indices for sample replicates and OTUs.

summarized_selectivity <- selectivity_data %>% 
  
  # Mean Si for all OTUs for sample replicates.
  dplyr::group_by(NODE_NAME, Class, Order, Family, Genus, Species, OTU, SAMPLE_date) %>% 
  dplyr::summarise(Selectivity_index= mean(Selectivity_index)) %>%
  
  # Merge OTUs at suitable level. (Here, Class)
  dplyr::group_by(NODE_NAME, Order, SAMPLE_date) %>% 
  dplyr::summarise(Selectivity_index= sum(Selectivity_index)) %>%
  
  #Filter uniportant contributions
  dplyr::filter(Selectivity_index > 0.00) %>% 
  dplyr::group_by(NODE_NAME, SAMPLE_date) %>% 
  dplyr::summarize(Selectivity_index=Selectivity_index/sum(Selectivity_index), Order)

head(summarized_selectivity)
rm(environment, paired_data, zoopl_16S)

1.2 Biomasses

The population biomasses of zooplankton and phytoplankton were retrieved from the Swedish national pelagic phytoplankton and zooplankton monitoring with sampling intensity ranging between monthly and weekly samples(SMHI 2018). Individual body masses of phytoplankton and zooplankton were retrieved from COMBINE guidelines for plankton monitoring in the Baltic Sea(HELCOM 2017b). We calculated daily biomass estimates over one year by linearly interpolating data from samples taken between 2007 and 2018.

Zooplankton <- readRDS("Data/Zooplankton_shark.rds")
Phytoplankton <- readRDS("./Data/Phytoplankton_shark.rds")

Daily interpolation of biomass

We calculated daily biomass estimates by linearly interpolating data from samples taken between 2007 and 2018. First we define the function for daily interpolation:

dailyInterpretation <- function(data, taxa="Genus") {
    
  require(tidyverse, xts)
  
  allDates <- seq.Date(
    #Creating a sequence of days for dataset
    min(as.Date(data$SDATE)),
    max(as.Date(data$SDATE)),
    "day")

  
  daydata <- data %>%
    
    # Turn into "matrix" format
    spread(key={{taxa}}, value=Value) %>%                 
    
    # Replaceing Nas with 0 for actual samples
    replace(is.na(.), 0) %>%                              
    
    # Adding the sequence of days to existing dataset
    full_join(x=data.frame(SDATE=allDates), y=.) %>%
    
    # Order acc to date
    arrange(SDATE) %>%
    
    # Linnear na approximation of all numeric variables in dataset
    mutate_if(is.numeric, funs(na.approx(.))) %>%                 
    
    # Back to long format
    gather(key = "Taxa", colnames(.)[-1], value = "Value") %>%    
    
    # Summarize days of each year
    mutate(Abundance = as.numeric(as.character(Value)),           
           DOY = as.numeric(strftime(SDATE, format = "%j")),
           Year = as.numeric(strftime(SDATE, format = "%Y"))) %>% 
    group_by(DOY, Taxa) %>% 
    summarize(Abundance = mean(Abundance))
  
  return(daydata)
}

Daily interpolation of zooplankton biomass

dip_zooplankton <- function(x) {

zoopl <- Zooplankton %>%
  dplyr::filter(Station == x,
        Depth %in% c(30,60),
         Year > 2006,
         Parameter == "Wet weight/volume",
         dev_stage_code != "NP",
         Genus %in% c(
           "Evadne",
           "Bosmina",
           "Acartia",
           "Temora",
           "Centropages",
           "Keratella",
           "Synchaeta",
           "Pseudocalanus"
           )) %>%
  group_by(SDATE, Genus) %>%
  
  summarize(Value=sum(Value)*30/1000) %>% # *30 To get the full water column, /1000 to get from mg to g.
  
  dailyInterpretation(taxa="Genus") %>% 
  mutate(NODE_NAME=Taxa, 
         Station=x)

  return(zoopl)
}

zooplankton_dip <- dip_zooplankton("BY31 LANDSORTSDJ") %>% 
  bind_rows(dip_zooplankton("BY15 GOTLANDSDJ")) %>% 
  bind_rows(dip_zooplankton("BY5 BORNHOLMSDJ"))
## Warning: `funs()` was deprecated in dplyr 0.8.0.
## ℹ Please use a list of either functions or lambdas:
## 
## # Simple named list: list(mean = mean, median = median)
## 
## # Auto named with `tibble::lst()`: tibble::lst(mean, median)
## 
## # Using lambdas list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
head(zooplankton_dip)

Daily interpolation of phytoplankton

dip_phytoplankton <- function(x) {


phytopl <- Phytoplankton %>% 
  
  filter(Station == x,
         Year > 2006,
         unit ==    "ugC/l") %>% 
  
  group_by(SDATE, Order, Depth) %>%
  summarize(Value=sum(Value)) %>%
  
  group_by(SDATE, Order) %>% 
  summarize(Value=mean(Value)*0.12) %>%
  # *0.12 to correct for depth *30, carbon content*4, and parameter *0.001 ug/l -> g/m
  
  dailyInterpretation(taxa="Order") %>%
  mutate(NODE_NAME=Taxa, 
         Station=x)
}

phytoplankton_dip <- dip_phytoplankton("BY31 LANDSORTSDJ") %>% 
  bind_rows(dip_phytoplankton("BY15 GOTLANDSDJ")) %>% 
  bind_rows(dip_phytoplankton("BY5 BORNHOLMSDJ"))
## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
## `.name_repair` is omitted as of tibble 2.0.0.
## ℹ Using compatibility `.name_repair`.
## ℹ The deprecated feature was likely used in the tidyr package.
##   Please report the issue at <]8;;https://github.com/tidyverse/tidyr/issueshttps://github.com/tidyverse/tidyr/issues]8;;>.
head(phytoplankton_dip)

1.4 Construct data with information of nodes

bodymasses <- tibble(NODE_NAME = c("Pseudocalanus",
                                   "Temora",
                                   "Centropages",
                                   "Acartia",
                                   "Evadne",
                                   "Bosmina",
                                   "Keratella",
                                   "Synchaeta"),
                     BODYMASS = c(1.69E-05,
                                  1.69E-05,
                                  1.69E-05, 
                                  9.50E-06,
                                  3.44E-05,
                                  3.44E-05,
                                  1.00E-06,
                                  6.00E-06))
nodenames <- Biomasses %>% 
  pull(NODE_NAME) %>% 
  unique() %>% tibble(NODE_NAME = .)

Nodes <- merge(bodymasses, nodenames, by=c("NODE_NAME"), all = T) %>%  
  mutate(BODYMASS = ifelse(NODE_NAME %in% bodymasses$NODE_NAME, as.numeric(BODYMASS), 1),
         ORG_TYPE = ifelse(NODE_NAME %in% bodymasses$NODE_NAME, "Animal", "Plant"))
         
head(Nodes)

1.5 Save all modified datasets

saveRDS(Links, "./Data/Mod/Links.rds")
saveRDS(Nodes, "./Data/Mod/Nodes.rds")
saveRDS(Biomasses, "./Data/Mod/Biomasses.rds")

2 Model Definition

2.1 Define FoodWebModel Function

The metabolic rate for each node was calculated based on metabolic scaling theory, where the metabolic rate Xi (J/s) is derived from species’ body masses Mi and an allometric scaling constant, i.e., a normalization constant a =17.17 for invertebrates and a = 18.47 for vertebrates. Bei represents the total population biomass per area unit (g/m2) and is multiplied by the metabolic rate (flux per gram biomass) to obtain a metabolic rate estimate at the population level. Metabolic rates are also adjusted to ambient temperature, where E is the activation energy, K is Boltzmann’s constant, and T is the absolute temperature in Kelvin. We used a bioenergetic model described by (Gauzens et al. 2019) to calculate energy fluxes (J/day/m2) between all nodes in the food web for each day of the year, and each of the three stations. The model builds on a steady-state assumption where the energy consumption of each species equals the energetic losses (Gi = Li). In this model, losses of each population (Li) were defined as the sum of metabolic losses for the population, and the energy flux lost to predation (Christensen and Pauly 1992). Gains for each population were calculated as the sum of all fluxes from prey to predator multiplied with the assimilation efficiency, i.e., conversion of consumed biomass into energy and an absolute prey preference constant. The efficiency was based on the functional group of the prey and put to 0.7755. The prey preference was based on the selectivity index calculated from DNA abundances scaled with the relative biomass abundance of the prey species.

#' Estimating fluxes of energy

#' For a given day and time, calculates fluxes in the food web, 
#' This is a wrapper function for fluxweb::fluxing, not intended for export.

#' @param nodes Nodes dataframe Should contain: BODYMASS, ORG_TYPE, and NODE_NAME
#' @param links DataFrame with links: Predator, Prey, Preference
#' @param biomasses Dataframe with biomass information: DOY, Abundance, NODE_NAME and Station
#' @param temperatures Dataframe with temperature: DOY and Temperature
#' @param doy Integer: What day of the year?
#' @param station Character string: Name of the station

#' @value data.frame

FoodWebModel <- function(nodes=Nodes,
                         links=Links,
                         biomasses=Biomasses,
                         temperatures=Temperatures,
                         doy=280,
                         station="BY31 LANDSORTSDJ") {
  
  
  # 1. Construct two support functions:
  
  #' Support function: Construct Adjacency matrix from edgelist

  #' @param links A dataframe with edges containing col1: Pred, Col2: Prey.
  #' @param weight Column name that contains whe edge weight.
  #' @value numeric.matrix
  
  edgelistToAdjacency <- function(links=Links, weight="Preference") {
    
    graph <- graph.data.frame(links)
    mat <- get.adjacency(graph, attr=weight) %>% as.matrix()
    
    adjacency <- mat %>%
      unlist() %>%
      as.numeric() %>%
      matrix(nrow=nrow(mat)) %>%
      t()
    return(adjacency)
    
    } # End of edgelistToAdjacency()
  
  
  #' Support function: Convert matrix to long format dataframe
  
  #' @param mat a matrix
  #' @value tibble. Col1: rownames(mat), col2: colnames(mat), col3: mat.
  
  matToDf <- function(mat) {
    
    DF <- list(Row = rownames(mat)[row(mat)] %||% row(mat),
               Col = colnames(mat)[col(mat)] %||% col(mat),
               value = mat) %>%
      map_dfc(as.vector) %>% 
      as_tibble()
    return(DF)
    
    } # End of matToDf()
  
  
  # 2. Calculate LOSSES and define EFFICIENCIES based on metabolic type and temperature
  temp <- temperatures %>% 
    dplyr::filter(DOY==doy) %>% 
    dplyr::pull(Temperature)
  
  boltz <- 0.00008617343  # Boltzmann constant
  tkonst <- 0.69/(boltz*(273.15+temp)) #Temperature metabolic constant
  
  nodes_tmp <-  nodes %>%
    dplyr::mutate(LOSSES = exp(-0.29 * log(BODYMASS) +
                                 dplyr::recode(ORG_TYPE ,
                                               "Animal" = 17.17,
                                               "ectotherm vertebrates" = 18.47,
                                               "Plant" = 1) -tkonst),

                  EFFICIENCIES = dplyr::recode(ORG_TYPE,
                                               "Plant" = 0.77,
                                               "Detritus" = 0.4,
                                               "Animal" = 0.906,
                                               "Bacteria" = 0.906))
  
  # Note: Detritus derived from:Gergs, R. and K. O. Rothhaupt (2008) in Freshwater
  # Biology 53: 2494-250.However, in practice all resources in this model are "Plants",
  # and all resources are "Animals", non ectothermal. All other numbers are not utilized. 
  
  
  # 3. Add BIOMASS to node data frame
  nodes_doy <- biomasses %>% 
    dplyr::filter(DOY==doy, Station==station) %>% 
    left_join(nodes_tmp, by="NODE_NAME")
 
  # 4. Rearrange dataframe to overlap with LINKS data.
  fluxnames <- links %>% 
    graph.data.frame() %>%
    get.adjacency(attr="Preference") %>%
    as.matrix() %>% 
    colnames()
  
  nodes_doy <- nodes_doy[match(fluxnames, nodes_doy$NODE_NAME),] %>% 
    as.data.frame()
  
  row.names(nodes_doy) <- nodes_doy$NODE_NAME
  
  
  # 5. Calculate energy fluxes usin the fluxweb packae
  
  flux <- links %>% 
    edgelistToAdjacency() %>% # Make square numeric matrix.
    fluxing(biomasses = nodes_doy$Abundance,
            losses = nodes_doy$LOSSES,
            efficiencies = nodes_doy$EFFICIENCIES,
            bioms.prefs = TRUE,
            ef.level = "prey",
            bioms.losses = TRUE)
  
  # 6. Clean up the output, and return it to dataframe formate.
  colnames(flux) <- fluxnames
  row.names(flux) <- fluxnames
  
  food_web <- flux %>%
    matToDf() %>%
    mutate(DOY=doy,
           Station=station,
           Predator = Col,
           Prey = Row,
           Flux = value) %>% 
    dplyr::select(Predator, Prey, Flux, DOY, Station)
  
  return(food_web)
  
} # End of FoodWebModel()

2.2 Computation and Calculations of metrics

Read in pre filtered data sets that the model is based on.

Links <- readRDS("./Data/Mod/Links.rds")
Nodes <- readRDS("./Data/Mod/Nodes.rds")
Biomasses <- readRDS("./Data/Mod/Biomasses.rds")
Temperatures <- readRDS("./Data/Temperatures.rds")

Calculate fluxes for for all days and stations

We calculated fluxes (J/day/m2) between all nodes in the food webs for each day of the year and each station. Here, the above defined FoodWebModel() function is applied to day 30:340, and appended to the three stations, one after the other.

Station <- c("BY31 LANDSORTSDJ", "BY15 GOTLANDSDJ", "BY5 BORNHOLMSDJ")
Doy <- 30:340
g <- expand.grid(Station, Doy)

full_model_df <-
  mapply(FoodWebModel,
         station = g[, 1],
         doy = g[, 2],
         MoreArgs = list(links=Links,
                         biomasses=Biomasses,
                         temperatures=Temperatures),
         SIMPLIFY = FALSE) %>%
  bind_rows() %>% 
  mutate(Flux=Flux*86.4) #Modify unit from j/m2/s to Kj/m2/day

Calculate metrics per day

For each node and day of the year, we calculated the total predation losses, normalized predation pressure, and total consumption. We also calculated annual food web metrics by summarizing all daily flux networks (J/year/m2).

Total

# Consumption per predator
Pred_flux <- full_model_df %>%
  mutate(NODE_NAME = Predator) %>% 
  group_by(NODE_NAME, DOY, Station) %>% 
  summarize(Consumption=sum(Flux))

head(Pred_flux)
Summary_statistics_stations <- full_model_df %>%
  mutate(NODE_NAME = Prey) %>% 
  
  # Calculate predation loss
  group_by(NODE_NAME, DOY, Station) %>% 
  summarize(Predation_loss=sum(Flux)) %>%
  left_join(Biomasses) %>%
  left_join(Pred_flux) %>% 
  mutate(Predation_pressure=Predation_loss/Abundance) %>% 
  gather(key = "Parameter", value = "Value",
         Abundance, Predation_loss, Predation_pressure, Consumption)

Summary_statistics_all <- Summary_statistics_stations %>% 
  group_by(NODE_NAME, DOY, Parameter) %>% 
  summarize(Value=mean(Value))

head(Summary_statistics_all)

2.3 Plotting

Define order and color scheme for plots

Order for food web plot

ord.pred.fw <- c("Temora", "Centropages", "Pseudocalanus", "Acartia",
              "Evadne", "Bosmina", "Synchaeta", "Keratella")

ord.prey.fw <- c("Synechococcales",
              "Chlorellales", "Pyramimonadales", # Chlorophytes
              "Prymnesiales", # Haptophyta, Hacrobia
              "Pyrenomonadales", #Cryptophyta, Hacrobia
              "Chromulinales", # Chrysophyceae, Ochrophyta, SAR
              "Thalassiosirales", "Chaetocerotales", #Baclillariophyta, Ochrophyta, SAR
              "Peridiniales", #Dinoflagellata, SAR
              "Nostocales")

color.pred.fw <- c("#8c510a", #Temora
                "#bf812d", #Centropages
                "#dfc27d", #Pseudocalanus
                "#f6e8c3", #Acartia
                "#c7eae5", #Evadne
                "#80cdc1", #Bosmina
                "#35978f", #Synchaeta
                "#01665e") #Keratella

color.prey.fw <- c("#3690c0", #Synechococales
                "#238443", #Chorellales 
                "#419753", #Pyramimonadales
                "#60AB64", #Prymnesiales
                "#7FBF75", #Pyrenomonadales
                "#9DD385", #Chromulinales
                "#FEA130", #Thalassiosirales
                "#FEB23F", #Chaetocerotales
                "#ef3b2c", #Peridiniales
                "#0570b0"  #Nostocales
)

Order for line diagrams

ord.pred <- c("Temora", "Centropages", "Pseudocalanus", "Acartia",
              "Evadne", "Bosmina", "Synchaeta", "Keratella")

ord.prey <- c("Synechococcales", "Nostocales",
              "Thalassiosirales", "Chaetocerotales", #Baclillariophyta, Ochrophyta, SAR
              "Peridiniales", #Dinoflagellata, SAR
              "Chlorellales", "Pyramimonadales", # Chlorophytes
              "Prymnesiales", # Haptophyta, Hacrobia
              "Pyrenomonadales", #Cryptophyta, Hacrobia
              "Chromulinales" # Chrysophyceae, Ochrophyta, SAR
              )

color.pred <- c("#8c510a", #Temora
                "#bf812d", #Centropages
                "#dfc27d", #Pseudocalanus
                "#f6e8c3", #Acartia
                "#c7eae5", #Evadne
                "#80cdc1", #Bosmina
                "#35978f", #Synchaeta
                "#01665e") #Keratella

color.prey <- c("#3690c0", #Synechococales
                "#0570b0", #Nostocales
                "#FEA130", #Thalassiosirales
                "#FEB23F", #Chaetocerotales
                "#ef3b2c", #Peridiniales
                "#238443", #Chorellales 
                "#419753", #Pyramimonadales
                "#60AB64", #Prymnesiales
                "#7FBF75", #Pyrenomonadales
                "#9DD385", #Chromulinales
                "#FEA130", #Thalassiosirales
                "#FEB23F", #Chaetocerotales
                "#ef3b2c") #Peridiniales


names(color.pred) <- ord.pred
names(color.prey) <- ord.prey


scale_color_pred <- scale_color_manual(breaks = ord.pred,
                                       values = color.pred)

scale_color_prey <- scale_color_manual(breaks = ord.prey,
                                       values = color.prey)

scale_fill_prey <- scale_fill_manual(breaks = ord.prey,
                                       values = color.prey)
scale_x_doymonth <- scale_x_continuous('Month',breaks=seq(0,365,by=30.5),
                                       limits=c(0,366),
                                       labels=c("Jan","Feb","Mar","Apr","May","Jun",
                                                "Jul","Aug","Sep","Oct","Nov","Dec"),
                                       expand = c(0,0))

Fig. 2 (inner) Annual Fluxes

fw_df <- full_model_df %>%
  
  # Calculate Sum fluxes of all days per year
  group_by(Predator, Prey, Station) %>% 
  summarize(Flux = sum(Flux)) %>%
  
  # Calculate Mean fluxes between station
  group_by(Predator, Prey) %>% 
  summarise(Flux = mean(Flux)) %>%
  
  # Make bipartite matrix format
  filter(Flux > 0)

fw_mat <- fw_df %>% 
  spread(Predator, Flux, fill=0) %>% 

  # Order the food web with defined orders
  relocate(ord.pred.fw) %>%   #Columns
  arrange(match(Prey, rev(ord.prey.fw))) %>%  #Rows
  column_to_rownames("Prey") %>%
  as.matrix()
## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
##   # Was:
##   data %>% select(ord.pred.fw)
## 
##   # Now:
##   data %>% select(all_of(ord.pred.fw))
## 
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
circos.clear()
chordDiagram(fw_mat, grid.col = c(rev(color.prey.fw), color.pred.fw),
             annotationTrack = c("grid", "axis"), directional = 1,
             direction.type = c("diffHeight", "arrows"), link.arr.type = "big.arrow")

Save plot:

Save summary matrix:

fw_mat %>% 
  as.data.frame() %>%
  rownames_to_column() %>% 
  write_csv("./Outputs/Tab.S1.csv")

Fig. 2 (outer) Fluxes per month

Prepare a food web matrix

fw_df <- full_model_df %>% 
  mutate(Month = month(as_date(DOY))) %>% 
  
  # Calculate Sum fluxes of all days per year
  group_by(Predator, Prey, Station, Month) %>% 
  summarize(Flux = sum(Flux)) %>%
  
  # Calculate Mean fluxes between station
  group_by(Predator, Prey, Month) %>% 
  summarise(Flux = mean(Flux)) %>%
  
  # Make bipartite matrix format
  filter(Flux > 0)

head(fw_df)

Filter and plot per month

Calculate monthly

fw_df %>% 
  group_by(Month) %>% 
  summarize(Flux=sum(Flux)) %>% 
  mutate(Fluxsqrt=sqrt(Flux))

Fig. 3 Selectivity indices

si_mat <- Links %>%
  
  group_by(Predator) %>% 
  summarize(Preference=Preference/sum(Preference),
            Prey = Prey) %>%
  mutate(Predator = factor(Predator, ord.pred),
         Prey = factor(Prey, ord.prey)) %>% 
  ggplot() +
  geom_bar(aes(x=Predator, y=Preference, fill = Prey), stat = "identity") +
  scale_fill_prey +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90))

si_mat

ggsave("./Outputs/fig.3.pdf",plot = si_mat, width = 5, height = 3)

Fig. 4 Plot Summary metrics

Fig. 4a Phytoplankton biomass

fig.4a <- Summary_statistics_all %>%
  filter(Parameter=="Abundance",
         NODE_NAME %in% ord.prey) %>% 
  ggplot(aes(DOY, Value)) +
  geom_line(aes(color=NODE_NAME)) +
  scale_x_doymonth +
  scale_color_prey +
  scale_y_continuous(expression(paste('Biomass(WW) (g / '~m^{2}~')'))) +
  labs(color='Prey') +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

fig.4a

Fig. 4b Zooplankton biomass

fig.4b <- Summary_statistics_all %>%
  filter(Parameter=="Abundance",
         NODE_NAME %in% ord.pred) %>% 
  ggplot(aes(DOY, Value)) +
  geom_line(aes(color=NODE_NAME)) +
  scale_x_doymonth +
  scale_color_pred +
  scale_y_continuous(expression(paste('Biomass(WW) (g / '~m^{2}~')'))) +
  labs(color='Predator') +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

fig.4b

Fig. 4c Predation losses (Secondary production)

fig.4c <- Summary_statistics_all %>%
  filter(Parameter=="Predation_loss",
         NODE_NAME %in% ord.prey) %>% 
  ggplot(aes(DOY, Value)) +
  geom_line(aes(color=NODE_NAME)) +
  
  scale_x_doymonth +
  scale_color_prey +
  scale_y_continuous(expression(paste('Secondary production (kJ / day /'~m^{2}~')'))) +
  labs(color='Prey') +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

fig.4c

Fig. 4d Predation Pressure

fig.4d <- Summary_statistics_all %>%
  filter(Parameter=="Predation_pressure",
         NODE_NAME %in% ord.prey) %>% 
  ggplot(aes(DOY, Value)) +
  geom_line(aes(color=NODE_NAME)) +
  
  scale_x_doymonth +
  scale_color_prey +
  scale_y_continuous(expression(paste('Predation pressure (kJ / day / g)'))) +
  labs(color='Prey') +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

fig.4d

Combine and save

fig.4 <- ggpubr::ggarrange(fig.4a, fig.4b, fig.4c, fig.4d)

ggsave("./Outputs/fig.4.pdf",plot = fig.4, width = 8, height = 5)

3 Model sensitivit to selectivity

To estimate the importance of zooplankton feeding preferences on annual primary, we resampled the selectivity index data for all predators using the “slice_sample” function in R base package (n=30, resampling with replacement). We also created a null model, assuming no prey preference, and no feeding niche differentiation between the predators, by randomly subsampling the selectivity dataset (n=30). The annual contribution of each primary producer to food web productivity (predation losses) was recalculated with the resampled selectivity indices and was compared with the null model, using linear regression of the log-transformed data.

3.1 Load data, summarize per order

summarized_per_order <- selectivity_data %>% 
  
  group_by(NODE_NAME, Order, SAMPLE_date, Sample.x) %>% 
  summarise(Selectivity_index = sum(Selectivity_index)) %>% 
  
  filter(Order %in% ord.prey.fw)
## `summarise()` has grouped output by 'NODE_NAME', 'Order', 'SAMPLE_date'. You
## can override using the `.groups` argument.

3.2 Subsample and bootstrap selectivity dataset

subsamplePreference <- 
  function(Subsample) {
    
    names_predator <-
      summarized_per_order %>% 
      spread(Order, Selectivity_index) %>%
      pull(NODE_NAME)
    
    names_prey <-
      summarized_per_order %>% 
      spread(NODE_NAME, Selectivity_index) %>%
      pull(Order)
  
    subsample <- 
      summarized_per_order %>% 
      spread(Order, Selectivity_index) %>% 
      group_by(NODE_NAME) %>%
      slice_sample(prop = 1, replace = TRUE) %>% 
      gather(key = Prey, value = Preference, Chaetocerotales:Thalassiosirales) %>%
      group_by(NODE_NAME, Prey) %>%
      summarise(Preference = mean(Preference)) %>% 
      mutate(Method = "grouped")

    random_predator <-
      summarized_per_order %>%
      spread(Order, Selectivity_index) %>%
      ungroup() %>% 
      slice_sample(prop = 1, replace = TRUE) %>% 
      mutate(NODE_NAME = names_predator) %>% 
      gather(key = Prey, value = Preference, Chaetocerotales:Thalassiosirales) %>%
      group_by(NODE_NAME, Prey) %>%
      summarise(Preference = mean(Preference)) %>% 
      mutate(Method = "random_predator")
    
    random_prey <- 
      summarized_per_order %>%
      spread(NODE_NAME, Selectivity_index) %>%
      ungroup() %>% 
      slice_sample(prop = 1, replace = TRUE) %>% 
      mutate(Order = names_prey) %>% 
      gather(key = NODE_NAME, value = Preference, Acartia:Temora, na.rm = TRUE) %>% 
      mutate(Prey = Order) %>% 
      group_by(NODE_NAME, Prey) %>%
      summarise(Preference = mean(Preference)) %>% 
      mutate(Method = "random_prey") 

    
    bind_rows(subsample, random_prey) %>% 
      mutate(Subsample = Subsample)
  }

3.4. Summarise and plot preference bootstrap data

Preference_summary <- 
  Links_I %>%
  group_by(Predator, Prey, Method) %>%
  dplyr::summarise(avg_Preference = mean(Preference),
                   uCI_Preference = Rmisc::CI(Preference)[1],
                   lCI_Preference = Rmisc::CI(Preference)[3])
## `summarise()` has grouped output by 'Predator', 'Prey'. You can override using
## the `.groups` argument.
Preference_plot <-
  Preference_summary%>%
  #filter(Method == "grouped") %>% 
  mutate(Predator = factor(Predator, ord.pred.fw),
         Prey = factor(Prey, ord.prey.fw)) %>% 
  ggplot(aes(fill=Prey, y=avg_Preference, x=Predator)) +
  geom_bar(stat = "identity", position = "dodge")+
  geom_errorbar(aes(ymin=lCI_Preference, ymax=uCI_Preference),position = "dodge", size=0.25) +
  facet_wrap(Method~.,
             labeller = as_labeller(c(`grouped` = "Model",
                                      `random_prey` = "Null-Model"))) +
  scale_fill_prey +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90)) +
  scale_y_continuous(expression(paste('Prey Selectivity Index')))

Preference_plot

3.5 Function: Estimate fluxes of energy

# Wrapper for FoodWebModel function

FoodWebModelI <- 
  function(nodes = Nodes,
           links = Links_I,
           biomasses = Biomasses,
           temperatures = Temperatures,
           doy=152,
           station="BY31 LANDSORTSDJ",
           iteration=1,
           method="grouped") {
    
  links_f <- links %>%
    filter(Subsample == iteration,
           Method == method)
  
  Model <- 
    FoodWebModel(nodes = nodes,
                 links = links_f,
                 biomasses = biomasses,
                 temperatures = temperatures,
                 doy = doy,
                 station = station) %>%
    mutate(Method = method, 
           Iteration = iteration) %>% 
    filter(Predator %in% ord.pred.fw,
           Prey %in% ord.prey.fw)
  
  return(Model)
}

Run the function for all stations, all days, stations, and iterations for the model and the null model.

Method <- c("grouped", "random_prey")   
Iteration <- 1:30
Station <- c("BY31 LANDSORTSDJ", "BY15 GOTLANDSDJ", "BY5 BORNHOLMSDJ")
Doy <- 30:340
g <- expand.grid(Method, Iteration, Station, Doy)

# Running this command may be slow, up to 1 h. Concider loading precalculated dataset.

full_model_iteration <-
  mapply(FoodWebModelI,
         method = g[, 1],
         iteration = g[, 2],
         station = g[, 3],
         doy = g[, 4],
         MoreArgs = list(links=Links_I,
                         biomasses=Biomasses,
                         temperatures=Temperatures),
         SIMPLIFY = FALSE) %>%
  bind_rows()

saveRDS(full_model_iteration, "./Outputs/full_model_iteration.rds")
full_model_iteration <- 
  readRDS("./Outputs/full_model_iteration.rds")

3.6. Summarize and plot output

summary_full <-
  full_model_iteration %>% 
  mutate(Flux=Flux*86.4) %>%  #Modify unit from j/m2/s to Kj/m2/day
  group_by(Iteration, Method, Station, Prey) %>% 
  summarise(Flux = sum(Flux)) %>% 
  group_by(Iteration, Method, Prey) %>% 
  summarise(Flux = mean(Flux)) %>% 
  group_by(Prey, Method) %>% 
  summarise(avg_Flux = mean(Flux),
            uCI_Flux = Rmisc::CI(Flux)[1],
            lCI_Flux = Rmisc::CI(Flux)[3])
## `summarise()` has grouped output by 'Iteration', 'Method', 'Station'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'Iteration', 'Method'. You can override
## using the `.groups` argument.
## `summarise()` has grouped output by 'Prey'. You can override using the
## `.groups` argument.
Flux_plot <-
  summary_full %>% 
  mutate(Prey = factor(Prey, ord.prey.fw)) %>% 
  ggplot(aes(fill=Prey, y=avg_Flux, x=Prey)) +
  geom_bar(stat = "identity", position = "dodge")+
  geom_errorbar(aes(ymin=lCI_Flux, ymax=uCI_Flux),position = "dodge", size = 0.25) +
  facet_wrap(Method~.,
             labeller = as_labeller(c(`grouped` = "Model",
                                      `random_prey` = "Null-Model"))) +
  scale_fill_prey +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90)) +
  scale_y_continuous(expression(paste('Energy Flux (kJ / year /'~m^{2}~')')))

Flux_plot

ggpubr::ggarrange(Preference_plot, Flux_plot, ncol = 1, common.legend = TRUE, legend = "right")

ggsave("./Outputs/Fig.S2.pdf")
## Saving 7 x 5 in image

3.7 Hypothesis tesing - What is the effect of selectivity?

summary_test <-
  full_model_iteration %>% 
  mutate(Flux=Flux*86.4) %>%  #Modify unit from j/m2/s to Kj/m2/day
  group_by(Iteration, Method, Station, Prey) %>% 
  summarise(Flux = sum(Flux)) %>% 
  group_by(Iteration, Method, Prey) %>% 
  summarise(Flux = mean(Flux)) 
## `summarise()` has grouped output by 'Iteration', 'Method', 'Station'. You can
## override using the `.groups` argument.
## `summarise()` has grouped output by 'Iteration', 'Method'. You can override
## using the `.groups` argument.
mod <- lm(log(Flux) ~ Prey*Method , data = summary_test)
Anova(mod)
plot(mod)

4 Session Info

sessionInfo()
## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur ... 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] car_3.1-0       carData_3.0-5   fluxweb_0.2.0   lubridate_1.8.0
##  [5] circlize_0.4.15 svglite_2.1.0   igraph_1.3.5    xts_0.12.2     
##  [9] zoo_1.8-11      forcats_0.5.2   stringr_1.4.1   dplyr_1.0.10   
## [13] purrr_0.3.5     readr_2.1.3     tidyr_1.2.1     tibble_3.1.8   
## [17] ggplot2_3.3.6   tidyverse_1.3.2
## 
## loaded via a namespace (and not attached):
##  [1] fs_1.5.2            bit64_4.0.5         httr_1.4.4         
##  [4] tools_4.2.1         backports_1.4.1     bslib_0.4.0        
##  [7] utf8_1.2.2          R6_2.5.1            DBI_1.1.3          
## [10] colorspace_2.0-3    withr_2.5.0         gridExtra_2.3      
## [13] tidyselect_1.2.0    bit_4.0.4           compiler_4.2.1     
## [16] cli_3.4.1           rvest_1.0.3         xml2_1.3.3         
## [19] labeling_0.4.2      sass_0.4.2          scales_1.2.1       
## [22] systemfonts_1.0.4   digest_0.6.29       rmarkdown_2.17     
## [25] pkgconfig_2.0.3     htmltools_0.5.3     dbplyr_2.2.1       
## [28] fastmap_1.1.0       highr_0.9           Rmisc_1.5.1        
## [31] rlang_1.0.6         GlobalOptions_0.1.2 readxl_1.4.1       
## [34] rstudioapi_0.14     shape_1.4.6         jquerylib_0.1.4    
## [37] generics_0.1.3      farver_2.1.1        jsonlite_1.8.2     
## [40] vroom_1.6.0         googlesheets4_1.0.1 magrittr_2.0.3     
## [43] Matrix_1.5-1        Rcpp_1.0.9          munsell_0.5.0      
## [46] fansi_1.0.3         abind_1.4-5         lifecycle_1.0.3    
## [49] stringi_1.7.8       yaml_2.3.5          plyr_1.8.7         
## [52] grid_4.2.1          parallel_4.2.1      crayon_1.5.2       
## [55] lattice_0.20-45     haven_2.5.1         cowplot_1.1.1      
## [58] hms_1.1.2           knitr_1.40          pillar_1.8.1       
## [61] ggpubr_0.4.0        ggsignif_0.6.4      reprex_2.0.2       
## [64] glue_1.6.2          evaluate_0.17       modelr_0.1.9       
## [67] vctrs_0.4.2         tzdb_0.3.0          cellranger_1.1.0   
## [70] gtable_0.3.1        assertthat_0.2.1    cachem_1.0.6       
## [73] xfun_0.33           broom_1.0.1         rstatix_0.7.0      
## [76] googledrive_2.0.0   gargle_1.2.1        ellipsis_0.3.2

5. References

Balvočiūtė, Monika, and Daniel H. Huson. 2017. “SILVA, RDP, Greengenes, NCBI and OTT How Do These Taxonomies Compare?” BMC Genomics 18 (S2): 114. https://doi.org/10.1186/s12864-017-3501-4.
Chesson, Jean. 1983. “The Estimation and Analysis of Preference and Its Relationship to Foraging Models.” Ecology 64 (5): 1297–1304.
Christensen, V., and D. Pauly. 1992. “ECOPATH II - a Software for Balancing Steady-State Ecosystem Models and Calculating Network Characteristics.” Ecological Modelling 61 (3-4): 169–85. https://doi.org/10.1016/0304-3800(92)90016-8.
Decelle, Johan, Sarah Romac, Rowena F. Stern, El Mahdi Bendif, Adriana Zingone, Stéphane Audic, Michael D. Guiry, et al. 2015. “PhytoREF: A Reference Database of the Plastidial 16S rRNA Gene of Photosynthetic Eukaryotes with Curated Taxonomy.” Molecular Ecology Resources 15 (6): 1435–45. https://doi.org/10.1111/1755-0998.12401.
Gauzens, Benoit, Andrew Barnes, Darren P. Giling, Jes Hines, Malte Jochum, Jonathan S. Lefcheck, Benjamin Rosenbaum, Shaopeng Wang, and Ulrich Brose. 2019. “Fluxweb: An r Package to Easily Estimate Energy Fluxes in Food Webs.” Edited by Sarah Goslee. Methods in Ecology and Evolution 10 (2): 270–79. https://doi.org/10.1111/2041-210X.13109.
HELCOM. 2017a. “Manual for Marine Monitoring in the Combine Programme of HELCOM.”
———. 2017b. “Manual for Marine Monitoring in the Combine Programme of HELCOM.”
Herlemann, Daniel P R, Matthias Labrenz, Klaus Jürgens, Stefan Bertilsson, Joanna J Waniek, and Anders F Andersson. 2011. “Transitions in Bacterial Communities Along the 2000km Salinity Gradient of the Baltic Sea.” The ISME Journal 5 (10): 1571–79. https://doi.org/10.1038/ismej.2011.41.
SMHI. 2018. “Swedish Database of National Monitoring.” http://www.sharkdata.se.